Onset of negative interspike interval correlations in adapting neurons 
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Negative serial correlations in single spike trains are an efTective method to reduce the variability 
of spike counts. One of the factors contributing to the development of negative correlations between 
successive interspike intervals is the presence of adaptation currents. In this work, based on a hidden 
Markov model and a proper statistical description of conditional responses, we obtain analytically 
these correlations in an adequate dynamical neuron model resembling adaptation. We derive the 
serial correlation coefficients for arbitrary lags, under a small adaptation scenario. In this case, 
the behavior of correlations is universal and depends on the first-order statistical description of an 
exponentially driven time-inhomogeneous stochastic process. 



I. INTRODUCTION 



Spike-frequency adaptation (SFA) is one of the main 
adaptation mechanisms in neural systems [U [2] . As its 
name implies, the effect defining SFA is the observation of 
a scaling in the input-output relationship between the in- 
jected current (or a stimulus property) and the firing rate 
of a spiking neuron, from an initial to a stationary map- 
ping [3:-6j. Several mechanisms can contribute to SFA 
(for example, depressing synapses [7 ); however, the most 
prominent mechanism accounting for SFA is the presence 
of (probably simultaneous) spike-related currents, which 
produce a negative feedback to the neuron in time scales 
ranging from tens to thousands of milliseconds [H [SHU] ■ 
Even when the full impact of these currents on neural 
coding is not completely understood, it is known that 
they contribute to the processing of static as well as tem- 
poral signals. For the processing of temporal signals, it is 
worth pointing out the high-pass filtering characteristics 
due to SFA [H [T2j [13] and its related sensitivity to input 
fluctuations [5] , the forward masking effect [TH [15] , the 
selectivity to complex stimuli [THl [IZ] , and the enhanced 
reliability of temporal coding |S] . 

For the case of static signals, given the absence of a 
temporal structure in the input, a neural system makes 
use of a rate code to represent them. A rate code is 
defined as the number of spikes in a certain tempo- 
ral window and it is completely described by the spike 
count statistics [l^- Spike- related adaptation currents 
have a twofold impact on this code. First, they modify 
the tuning curve between signals and responses or the 
input-output relationship mentioned above, which helps 
to match dynamic ranges [2HS]- Second, they introduce 
negative correlations between successive interspike inter- 
vals (ISIs) in stationary neural spike trains ITMT^ IT^ . 
While the first effect reflects the modification of the mean 
spike count, the second effect implies a strong change in 
its variance (and higher-order properties). This change 
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arises from the fact that the presence of negative cor- 
relations in a point process reduces the long-term vari- 
ability in the counting process that defines the rate code 
[20]. Taken together, both effects deeply affect the en- 
coding reliability [121 [HHH] • The presence of negative 
correlations also affects the coding capabilities of other 
related schemes; for example, coding of slowly varying 
signals through a mechanism called noise shaping |251 - 
\T7\ (demonstrated for negative correlations arising from 
a history-dependent threshold, but also valid for spike- 
related adaptation currents), or adaptation-based inde- 
pendent codes [28] . 

Correlated single spike trains constitute a nonrenewal 
point process. In neural systems, this kind of process is 
relatively ubiquitous [29H33] (for reviews, mostly based 
on negative correlations, see also [THl [IS])- In gen- 
eral, there is a coexistence of processes that evoke op- 
posite effects on the ISI correlations, and therefore on 
the counting statistics in single neurons (adaptation cur- 
rents and other regularizing processes such as synaptic 
depression and negative feedback versus filtering and in- 
put correlations, among others). Such interesting scenar- 
ios have attracted relative attention within the theoret- 
ical neuroscience community, and several studies focus 
on these statistics or related properties in different situ- 
ations [231 [241 [261 [281I33H37]. Related with our present 
work, we should note three approaches that have been de- 
rived recently: a population-based scheme for adapting 
ensembles [3H] (see also [2311M11211), a directed discrete 
representation for counting events with a general internal 
structure [35] (see also [IHIEI]), and a general framework 
for nonrenewal processes as a hidden Markov model [H] . 
In particular, our work can be framed within the general 
approach obtained by van Vreeswijk [42] . 

In this work, based on the results we have found in a 
previous study about the first-passage-time (FPT) prob- 
lem in an exponentially driven Wiener process [43 , we 
derive how the resulting negative correlations arise in the 
spike train of a dynamical, although simple, process re- 
sembling the basic features of a spike-related adaptation 
current added to a spiking neuron in the presence of 
fast additive fluctuations. The expressions we find are 
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strictly valid in a slight adaptation regime, where the 
complete FPT density is not necessary and perturbation 
techniques can be applied, so they remain valid for other 
dynamical models (with the spike-related adaptation cur- 
rent considered here, or similar). In this way, the onset of 
correlations due to adaptation is general across different 
models, provided the additive noise is fast. In the final 
part of the work, we use the statistics of the FPT problem 
and the emerging correlations due to adaptation to assess 
how the spike count variance decreases in comparison to 
an equivalent neuron without adaptation, in an asymp- 
totic limit. This reduction in the spike count variance 
underlies the improvement in the decoding performance, 
and we show how the dependence on the correlations and 
on the intrinsic variability reduction due to the inhomo- 
geneous driving shape the spike count variance reduction 
for different spiking frequencies. 

II. THEORETICAL FRAMEWORK 

A. Basic model 

We consider an integrate-and-fire (IF) neuron, where 
subthreshold dynamics of the membrane potential V is 
governed by 



which is widely expressed in neurons exhibiting SFA [x{t) 
would represent the calcium concentration in a current- 
based scheme]. Without mathematical loss, we can set 
9a — Cm/ Til and use a to control the strength of the 
adaptation current. 

Since /adapt(i) = -(Cm/ra) x{t), from Eq. ^ it follows 
that, between the (n)th and the (n + l)th spikes, the 
evolution of the adaptation current is given by 

exph(i-t„)/ra], (3) 

where £„ represents the state of the adaptation process 
X immediately after the spike time t^- 

Without a random component, the deterministic dy- 
namics of the membrane potential in the IF model with 
an adaptation current, Eqs. ([T]) and ([s]), evolves along 
a prescribed trajectory and no correlations emerge since 
each ISI is a replica of itself (however, even in the deter- 
ministic regime, perturbations propagate in the sequence 
of ISIs, which induces correlations '37'). We introduce 
a stochastic component in the system through external 
noise. In particular, we assume that the external current 
is given by a constant deterministic component and an 
additive Gaussian white noise representing fast external 
fiuctuations 
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External, /cxt, as well as internal currents, f{V) and 
^adapt, drive the membrane potential whenever V < Vthr- 
The internal current f{V) takes into account different 
interspike phenomena, such as leakage or spike initiation 
onset. The simplest models are the perfect [f{V) — 0] 
and the leaky [f{V) = —ghV] IF neurons, where only 
leakage is considered (the perfect IF model corresponds 
to no leakage). The subthreshold dynamics is supple- 
mented by a threshold condition, which simplifies the 
highly nonlinear process of a spike excursion: whenever 
the potential reaches Vthr a spike is defined and immedi- 
ately, the membrane potential is set to a reset condition 

Several types of adaptation currents have been char- 
acterized by experiments [8141 0|. According to previous 
theoretical studies [H [13 [H] , we model the adaptation 
current as a process x{t) that decays during spikes and is 
incremented when a spike event occurs. In particular, we 
consider the adaptation current as /adapt (0 — ~9a. x{t), 
where the interspike dynamics for the adaptation process 
is given by 



dx 
'dt 



X 



(2) 



and a fixed increase a > is evoked at all spike times 
^sp^ 2;(isp) — >■ x(tsp) -fa. This model could be considered 
as an idealization of the Ca^+-dependent current. 



(4) 



where (C(t)) = and (^(i)C(i')) = 2D5{t' --t). Therefore, 
the subthreshold dynamics between the (n)th and the 
(n + l)th spikes is determined by 

^ = g{V) + M - - cxp[-(t - i„)/ra] + ^(t). (5) 

dt Ta 

Different (IF) models will have different g{V) func- 
tions. In particular, g{V) = for the perfect IF model 
and g{V) ~ —V/t^, with = C^jgL, for the leaky 
IF model. The threshold condition (at spike time tn+i) 
resets the membrane potential to K and updates the 
adaptation state to En+i — £n exp(— Tn/ra) -I- a, where 
Tn = ^n+1 — IS the (n)th ISI. 

In Fig. [ija) we show a typical realization of the mem- 
brane potential, V{t), and the adaptation process, x{t)^ 
for the system described above. As shown in Fig. [l](b), 
this model exhibits SFA, where a step input (bottom) 
induces an initial firing rate /q which decays to a lower 
steady-state firing rate fss (top), with some typical time 
scale. Spike times tn defining the spike train, Tsp(i) — 
'^5{t — in) [see top of Fig. [l](a)], also establish the se- 
quence of FPT processes, {. . . , t„_i, Tn, t„+i, . . . } = {t^} 
[Fig. [T]jc)]. The ongoing (initial) state of the adaptation 
process, £„, is a history-dependent random variable turn- 
ing the sequence {rn} into non-Markovian. 

However, the bidimensional state s = (e, r) , defined 
at spike times, supports a Markovian process since the 
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FIG. 1; (Color online) Neuron model with adaptation current, (a) A system realization, V{t) and x{t), and the resulting spike 
train Tsp{t). (b) The system develops spike- frequency adaptation (top), usually characterized by the temporal evolution of the 
firing frequency, f{t), stimulated with step currents (bottom), (c) A sample trace of the adaptation process x{t) and the basic 
elements used to describe the stochastic nature of the system, e and r. (d) The hidden Markov model, composed by e and the 
observable t, used for the analysis. 



state at (spike) time is completely characterized from 
the knovifledge of the state at (spike) time tn-i- given En-i 
and Tn_i, £n is defined by the deterministic relationship 
En = £n-i exp(— Tn-i/Ta) + «, and Tn is given by the FPT 
problem of a certain stochastic process with an exponen- 
tial time-dependent drift ( — en/'T'a) exp[— (t — tn)/Ta]. In 

particular, for the perfect (leaky) IF model, the underly- 
ing stochastic process is a Wiener (Ornstein-Uhlenbeck) 
diffusion process. In mathematical terms, the transition 
probability density is given by 

/(s„|Sn-l,Sn-2, • • • ) 
= /(^n, Tnlen-l, Tn-l) 

= (5{e„ - [£„_i exp(-Tn-l/Ta) + Q;]} 0(rn|£n), (6) 

where </>(Tn|e,i) is the FPT density function associated to 
the (n)th ISI, which depends exclusively on En and not 
on previous outcomes (tn in the above notation for the 
drift just set the initial time). This Markovian process is 
represented in Fig.jljd) and constitutes a hidden Markov 
model. Based on Eq. (|6]), the key elements necessary to 
analyze this stochastic system are the statistics of s and 
the time-inhomogeneous FPT density function. 



B. Statistics of the (initial) adaptation strength 

For one-dimensional systems as the IF models, the so- 
lution to the FPT problem is given as an expansion in 
terms of the strength of the exponential drift [43, ,44] 



0(r|e) = 5]e'</)Kr). (7) 

i=0 

Given the transition density, Eq. ([6]), the equilibrium 
probability density for e, Pcq(e), should satisfy [12] 



Poq(en) = J /e(£nkn-l) /3cq(£n-l) dSn-l, (8) 

where the transition density between substates is 



n 7 '^n l^n-1 1 "^n- 

(9) 

Replacing Eq. Q and the e-expansion for 0(Tn-i |£n-i), 
Eq. ([7|, in Eq. ([9f we obtain a self-consistent integral 
equation, hard to tackle analytically. However, from this 
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Laplace transform of the jth term in the expansion for 
the FPT solution. The previous relationship relates 
unconditional moments (subindexes are irrelevant), and 
therefore it represents a set of infinite algebraic equations 
for the (infinite) moments (e™). 

In the slight adaptation regime, a ^ but small, it is 
easy to see that the moments (e™) ~ 0{a"^). In partic- 
ular, the first two moments read 



£ = 



a 



[l-0O^(l/^a)]' 

[l + 0o^(l/ra)] 
[l-0^(l/ra)] [l-^^(2/ra)] 



(11) 
(12) 



The a-normalized properties, (1/Q;)(e) and {\/a'^){{e— 
(e))^), are shown in Fig. l2[a) in black and gray arrows, 
respectively (right margirij^ for the case of the perfect IF 
neuron model. As expected, they agree with the results 
obtained from simulations in the asymptotic limit (sym- 
bols plus error bars: mean plus standard deviation of £„; 
colored histogram: logarithm of the density). 

For small a values, there is an additional interesting 
viewpoint to derive (e) , also valid for any neuron model. 
Multitrial experiments usually start from rest, where we 
assume that there is no adaptation, Eq = 0- Conditional 
to this fact, the (n)th event satisfies 



FIG. 2: (Color online) Statistics of e for the perfect IF model, 
(a) Density of normalized e as a function of the spike number 
(color bar in logarithmic scale) . Symbols and error bars indi- 
cate mean and standard deviation, respectively. The analyti- 
cal expression for the mean as a function of the spike number, 
Eq. ( 16 1, is represented by the dotted (blue) line. The analyt- 



ical expressions fo r th e stati ona ry properties (mean/standard 
deviation), Eqs. (Ill and (12 1, are indicated by arrows at 



the right (black/gray). Numerical results were obtained for 
II = 0.10 ms-\ D = 0.05 ms-\ Khr-'K = 1-0, = 100.0 ms, 
and a = 0.01 (number of trials, iVtriais = 10^)- (b) The nor- 
malized stationary mean as a function of a (semilogarithmic 
plot). The analytical result is represented by a continuous 
(red) line. Nonlinear effects for large a are fitted by a (dashed- 
dotted) line, indicating that linear regime remains valid up to 
a ~ 0.1. Inset: plot of the (not normalized) stationary mean 
as a function of a. Same parameters as in (a), except that 
Nt-Azis = 10^ for each point. 



integral equation it is relatively simple to find a relation- 
ship between the moments, which reads 



(^^n ) = j Pcq(£n) ds-a 

= T E'^f(i/^a) (41)(10) 

i=l ^ j=0 

where (^) is the binomial coefficient and 4>f{s) is the 



e^ = a[l + Y,l[(^M-r^/re.)]- (13) 

i=i j=i 

The first moment is given by 



(e„) =a{l + Y^iU cxp(~rj/T,))]. (14) 
i=i j=i 

For a weak adaptation process we consider that (£„) 
0{a), and then the conditional probabilities required in 
the above expression are well approximated by their zero- 
order expansion, friTj^/Tn-i, Tn-2, • • • ) — 4'oiTn) for all n. 
In this limit it is easy to see that 

n-l n-l 

(en) = a [1-f ^ J|(exp(-rj/ra))0j 
i=i j=i 

= a $][0o^(l/ra)]\ (15) 

where (•)0; indicates the average with respect to the func- 
tion 0i(r). This geometric series is readily obtained: 



(£„) = a 



1 - [0o^(l/ra)]" 

[l-0^(l/^a)] ' 



(16) 
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and converges asymptotically to a/[\ — 0q (l/ra)], when- 
ever 0Q(l/ra) < 1, in concordance with the previous 
analysis. The exponential growth of the normalized 
adaptation strength (!/«)(£„), predicted by Eq. (16), is 
shown in Fig. [2ja) as a function of the spike number 
n (blue dotted line), together with the results obtained 
from numerical simulations (symbols). The agreement is 
remarkable for this case (simulation results were obtained 
with a = 0.01). 

In order to determine the range of a where the approx- 
imation remains valid, we run several simulations and 
calculate the asymptotic (equilibrium) (e) for different a 
values. In Fig. [2|^b) we show the normalized asymptotic 
mean value as a function of a (inset: not normalized 
mean value). The average exhibits the linear behavior 
indicated by the preceding results up to a ~ 0.1, which 
represents an adaptation of about 10% [(/o — fss)/ fo, see 
Fig-Sc)]. 



III. RESULTS 



where dsi represents dsidri and X'si its integration do- 
main. 

At lag k = 1, the transition probability density be- 
tween states is given directly by Eq. (6|. The onset of cor- 
relations is characterized by the smallest order in a which 
produce finite SCC values. This order coincides with 
the small adaptation regime, and therefore, the linear 
e-expansion suffices for relevant expressions in Eq 
(TnTn+i) ^ 0{a). Consequently, we obtain 



(191, 



[TnTn+l) 



l<t>o 
M 

a 



ds 



Ji/r. 



,(20) 



A. Onset of correlations 



where indicates the average with respect to the func- 
tion (piir). In order to keep the linear order in the nor- 
malization required for the SCC, we have (r)^ ^ 0{a) 
and (Ar2) - 0(0): 



Correlations in nonrenewal point processes are usually 
characterized by the serial correlation coefficient (SCC), 
which is defined by 



:Ar2)(Ar„V) 



(17) 



where k indicates the lag between successive ISIs, brack- 
ets indicate ensemble average, and Ar^ — (ti — (ti))^ 
is the variance. Once the system reaches the station- 
ary conditions (the firing frequency is adapted) , the SCC 
simplifies to 



(r)2 = (r)^^+2(£)(r)^„(r)^,, (21) 
(Ar^) = (Ar2)^„. (22) 



In Eqs. (20) and (21), the evaluation of (e) should be 
performed in the linear regime. Taking into account the 
results obtained in the previous section, {s)/a = 1/[1 — 
</>^(l/Ta)], the first SCC reads 



(Ar2 



(18) 



For the hidden Markov model defined by Eq. (|6|, 
(TnTn+k) reads 



Pi 



[l-0o^(l/ra)] (Ar2)^„ 



0o(l/^a)(T)0o + 



ds 



Jl/Ta 



(23) 



TnTn+k) = / Tn^'n+k /(^n+k , ^n+k kn , Tn) 

X 0(Tnkn) Pcq(£n) dSn dSn+k,(19) 



To evaluate the SCC at higher lags, k > 1, we need 
to obtain the transition probability density between the 
states n and n + k which, based on the hidden Markov 
model, reads 



/(£n-|-k, Tn-l-klen, Tn) — 



£n+k, Tn+k, • ■ • , £n+l, Tn+l\Sn, Tn) 



ds^^ ]^/(en+i,Tn+i|en+i-l,T„+i-l), 



(24) 
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FIG. 3; (Color online) Normalized serial correlation coefficients (SCC), pk = pk/ct, for the perfect IF model, (a) Theoretical 
expression for the normalized first SCC as a function of driving parameters, p, and D (ms~^) [see Eq. ( 23 1]. (b) Numerical results 
(symbols) for the normalized first SCC as a function of p (ms~^) and different D values (ms"'^). Corresponding theoretical 
expressions are indicated by continuous lines, (c) Factor (^Q(l/Ta), which defines the relationship between successive SCC 
values, as a function of the driv ing parameters, (d) Theoretical prediction of normalized SCC at higher lags as a function of 
driving parameters [see Eq. (26l]. Continuous, dashed, and dashed-dotted (colored) lines represent pi, p2, and pz, respectively, 
at different noise intensities (black, red, blue, and orange: D = 10°, 10"^ 10"*, and 10"*^ ms , respectively). Lines are 
embedded in their respective continuous (semitransparent) surfaces, (e) Numerical results for the normalized SCC at higher 
lags, pi, p3, and ps as a function of p (ms~^, D — 0.05 ms~^). Continuous lines represent the corresponding theoretical 
expressions. Dashed black line indicates 4>q{1/ts) (right scale), (f) Numerical results for the normalized SCC as a function of 
the lag for the four cases indicated with arrows in (e). Dashed lines indicate theoretical expressions. Ratios between successive 
SCCs are given by the value of <^o(l/''a) at each point [dashed line in (e) at the different values of p\. Numerical results were 
obtained from A''isi = 10^ consecutive ISIs, error bars were estimated from Nrepet = 10 repetitions. Remaining parameters: 
Vthr - K = 1.0, Ta = 100.0 ms, and a = 0.1. 



where we have simplified the notation to ds^'^"^-' = 
dsn_|_i . . . dsn+k-i and 'Ds^^^^'> to the corresponding in- 
tegration domain. Each of the factors under the product 



symbol has the structure given by Eq. ([6]). After some 
calculus [it is convenient to leave unevaluated all inte- 
grals in Tj in the transition density from to Sn+k, for 



7 



posterior evaluation in Eq. (19)], the average between 



successive ISIs is given, up to first order in a (and/or 
by 



(a) 



o.oa 




nT^n+k/ 



k-1 



(l + ^)(r)^„ + (r)^„^[</.^(l/ra)]'^-' 



/-r Mk-l ^'/'0(-'^) I 
— [0o(l/ra)] ^^Jl/r„ 



.(25) 



Replacing the stationary value for (e), and after normal- 
ization, it is relatively simple to prove that SCC values 
at superior lags are given by 



ik-l 



Pi 



k> 1. 



(26) 



For the perfect IF model, the lowest orders in the ex- 
pansion corresponding to the solution of the FPT prob- 
lem, Eq. ([T]), are given [53], and therefore, the expressions 
for pk can be explicitly evaluated. 

In Fig. |3ja) we show the normalized first SCC, pi — 
Pi/a, as a function of the parameters governing the dy- 
namics of the perfect IF neuron model, fi and D. In 
Fig. [3](b) we show the theoretical expression for the nor- 
malized SCC pi as a function of the driving force /i, for 
different noise intensities, and compare it with numeri- 
cal results. The agreement between both curves is re- 
markable, since the limit of small adaptation is satisfied 
(numerical results were obtained with a = 0.1 and then 
normalized) . 



As expressed by Eq. (|26|, the SCC at higher lags, pk for 

k = 00 (l/'^a) Pk-1- 

c)], SCC at higher 



k > 1, have a geometric structure: 
Since < (^o (1/T"a) < 1 [ 

see Fig. 

lags are scaled versions of the first SCC. In Fig. Isjjd) we 
show the first three normalized SCC, pk — Pk/ce, as a 
function of the parameters p. and D. As stated by the 
geometric relationship, the exact scaling between them is 
given according to the precise location in the parameters 
space [equivalent point in Fig. [sjjc)]. The scaling can be 
better appreciated in Fig. [sjje), which simply represents 
a section of Fig. [sjjd) along a particular D. In this 
case, theoretical (solid lines) as well as numerical results 
(symbols) for pi , , and ps are represented as a function 
of the driving force p (p2 and p4 are omitted for the sake 
of clarity). The dashed line shows the respective section 
of Fig. (scale in the right margin), which governs 
the scaling between consecutive SCCs. Actually, the 
geometric structure represents an exponential decay of 
the SCC as a function of the lag. This can be observed 
in Fig. [sjjf) for different p values selected in Fig. [sjje). 
The exponential decay is given from the scaling factor 
(/)q (1/Ta), obtained from the intersection of the dashed 
line in Fig. [sf^e) and the particular value of p considered. 
For example, p2 and p4 were selected so pi were 
approximately the same for both cases. However, the 
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FIG. 4; (Color online) Normalized serial correlation coeffi- 
cients (SCCs), pk = pk/ct, for the leaky IF model[(/(V) = 
^V/Tm]- (a) and (b) analogous to Figs. [Sje) and ^i). Pa- 
rameters: Tin = 10.0 ms, D — 0.01 ms~^, Vthr — VI- = 1.0, and 
Ta = 100.0 ms. 



scaling factor (^Q(l/ra) is higher for p4 than for p2, so 
the decay is correspondingly slower. 

The onset of correlations, as characterized by the 
SCC, has a general structure, Eqs. (23) and (26), that 
relies critically on two factors: the Laplace transform 



of the unperturbed distribution, (/)q (s), and the linear 
correction to the mean introduced by the exponential 
temporal drift, (t)^^. As shown in |43| . for small 
intensities the main effect of the exponential drift on the 
FPT statistics of a perfect IF model is to change the 
mean consistently to what we have found in this work. 
For other IF models, the complete FPT problem with 
an exponential temporal drift can be addressed with a 
similar procedure |44j . However, even when appealing 
to set out the formalism, the complete statistics is not 
required for computing the SCC and so we can proceed 
with simpler approaches. For example, to compute 
the linear correction to the mean FPT due to the 
exponential temporal drift in generic IF models, we can 
use the results obtained by Lindner using a perturbation 
scheme [45l |46]. To illustrate the generality of the 
results we have found, in Fig. [4] w e show the normalized 
SCC, pk, predicted by Eqs. (|23|) and ((26|), using the 



analytical results, (f>Q{s) and (r)^^ , obtained by Lindner 
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FIG. 5: (Color online) Loss of linearity, (a) First SCC, pi, 
as a function of the driving parameter p in a perfect IF neu- 
ron model, for different adaptation strengths, a. In all cases, 
D = 0.02 ms"\ Khr -Vr = 1.0, and Ta = 100.0 ms. (b) As 
the value of a increases, the stationary firing rate decreases 
[asymptotic value of f{t) at the right margin] and correspond- 
ingly, the adaptation effect is more prominent. In all cases, 
the driving drift is /i = 0.10 ms~^ (which corresponds to a 
firing rate of 100 Hz in the absence of adaptation, a — 0), 
applied as a step function at time t = [see Fig. [ijb)] . Re- 
maining parameters as in (a). 



for the leaky IF neuron model [l^, and compare them 
to numerical simulations. As expected, theoretical and 
numerical results agree. 



in Fig. [5ja) as a function of ^ or D] , but only a rough 
estimate of the correct values of the correlations (cyan 
symbols and Hne in Fig. |5| . 



C. Spike-count variance reduction 

The presence of negative correlations affects the spike- 
count variance. To analyze this effect, it is convenient to 
introduce the Fano factor, FFt, which relates the mean 
and the variance of the spike counts observed in a tem- 
poral window of length T, 
as the ratio 



{nx) and (Ari^) , respectively. 



FF, 



(AnI) 
(tit) 



(27) 



For point processes, the asymptotic behavior of the 
Fano factor reads [201 



FF^ = hmr-^oo FF^ = CV^ 1 



OO > 

k=i / 



(28) 



where CV is the coefficient of variation, defined as the 
ratio between the standard deviation and the mean of 
the unconditional ISI statistics. 



CV: 



(29) 



B. Loss of linearity 

In the previous section we have shown that the on- 
set of correlations has a specific structure and scales lin- 
early with the adaptation strength, Eqs. (23) and (26 1 



This scaling enables us to consider normalized SCC, pk, 
as shown in the preceding figures. However, for a large 
enough value of a, higher-order effects become important 
and these equations are no longer applicable. In partic- 
ular, given that the expression for pi, Eq. (231, does not 



saturate we expect that higher-order effects oppose the 
linear growth. In Fig. [sjja) we show the (not normal- 
ized) first SCC, pi, for different adaptation strengths. 
As expected, for small values of a the theoretical ex- 
pression, Eq. ([23]), properly accounts for the linear scal- 
ing. Specifically, correlations grow linearly up to approxi- 
mately a ~ 0.25 [red symbols and line in Fig.[5ja)], which 
represents a frequency adaptation of about 20% [see red 
line in Fig. [sjjb)]. This limit is slightly better than that 
expected from the upper limit of the linear scaling in the 
stationary mean adaptation strength, (e) [see Fig. [2]jb)]. 
For larger values, higher effects are non- negligible and 
affect the correlations in the predicted manner. For a 
realistic SFA (adaptation of about 50%), the analyti- 
cal expressions provide a qualitative agreement regarding 
the dependence on the parameters [shape of the curve 



Combining the preceding equations and given that 
(ut) — T/{t), the spike-count variance reads, in the 
asymptotic limit. 



(An^) = 



(At2 




(30) 



The linear growth in T of the spike-count variance is 
a characteristic of a diffusive process (and the reason 
for the usefulness of the Fano factor). Inasmuch as the 
asymptotic limit is reached, it is useful to analyze the 
preceding factor, which we denote (An^), 



(Ar2 



2^Pk 



(31) 



k=l 



Equation (31) highlights two contributions to the 



spike-count variance: a contribution from the FPT statis- 
tics related to a single spiking process, (At^)/(t)^, and 
a contribution from the ISI correlations provided by the 
entire spike train, (1 -I- 2^5^^^ pk)- 

The effect of the adaptation strength a on the spike- 
count variance is twofold; it changes the ISI statistics as 
well as the correlations between the ISIs. In the slight 
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where (e) is given by Eq. (11). Replacing the com- 



plete expressions for both factors [Eqs. (32 1 and (33)] in 



Eq. ( 31 ), and keeping the first order in a [{An'^) 
we obtain 



FIG. 6: (Color online) Spike-count variability in the perfect 
IF neuron model, as a function of the firing frequency. The 
neuron without adaptation (a — 0) shows a spike-count vari- 
ance that does not depend on the firing frequency [theoretical 
(T) and numerical (A'^) results are represented by a black thick 
line and open squares, respectively]. For a ^ {a = 0.05 and 
a = 0.10 in red and green, respectively) there is a reduction 
in the variance, which also is independent of the firing rate 
[theoretical (T) and numerical (A*') results are shown by thick 
lines and empty symbols, respectively]. By shuffling the orig- 
inal spike train, from which the spike counts are drawn, we 
destroy the correlations. For the shufHed spike train, theoret- 
ical (r sh.) and numerical (A'^ sh.) results are represented by 
hidden lines and semifiUed symbols, respectively. Parameters: 
Vthr - Vr = 1.0, Ta = 100.0 ms, D = 0.02 ms-\ The driving 
parameter fj, was varied from 0.005 to 0.2 ms~^. Numerical 
results are presented as a function of the firing frequency ob- 
served from the mean spike counts. Theoretical results are 
shown as a function of the firing frequency computed from 
= (r)^f, -|- (e)(T)</,i. Numerical procedure: to compute 
the spike-count statistics, the length of the temporal window 
T was varied for different values of jj., checking in all cases 
that the asymptotic regime was reached. Numerical results 
were obtained from a single spike train with a temporal length 

Ttotal = T A'^aample, with iVsamplo = 100. The ValuC of A^samplo 

is low, due to the computational cost of the shuffling pro- 
cedure. To improve the accuracy, results are presented as 
averages over A^rcpot = 500 repetitions. 



adaptation regime we consider, the contribution to the 
spike-count variance due to the correlations is a linear 
term in a. Explicitly, it is easy to show that 



k=l 



Pi (a) 



(32) 



where pi{a) is given by Eq. (23). To be consistent with 



this description, the term provided by the ISl statistics 
should be linearized, and reads 



{r)l [l-</'o^(l/Ta)](r)3„ 



/0i 



(r>0o 1-'^0^(1/Ta) 

X ((r)^„ - 3(r)^„0o^(l/ra) - 2^^Ji/..) 



.(34) 



It is clear that the independent term on the right-hand 
side of Eq. (34) corresponds to the asymptotic spike- 



count variance (normalized by T) of the neuron without 
adaptation {a = 0). Even when general in the regime 
we consider, the behavior of this rather complicated ex- 
pression for a 7^ is not obvious at all. The expression 
simplifies enormously for the case of the perfect IF neu- 
ron model. For this neuron, we have 



(An|) 



2D 



AD 



(T4hr - V,Y (Vthr - V,Y- 



a. 



(35) 



Surprisingly, this expression for the spike-count vari- 
ance reduction does not depend on the driving param- 
eter ^, which sets the firing frequency. In Fig. |6] we 
show this behavior for different values of the adaptation 
strength. The case a = corresponds to a perfect IF neu- 
ron model without adaptation (black line and squares). 
As a increases, the spike-count variance decreases, and 
this effect does not depend on the firing frequency set 
by varying /i. For a = 0.05, theoretical and numerical 
results agree (red thick lines and empty circles, respec- 
tively); whereas for a = 0.10 the flat reduction is numer- 
ically observed, but the appropriate variance decrease is 
just approximated by the theoretical prediction (green 
thick lines and empty triangles). The same holds true for 
higher adaptation strengths (not shown). As expected, 
since the expression for pi [Eq. (|23[ )] and the exponential 
structure for higher lags [Eq. (|26[)] are approximate, the 
infinite sum in Eq. ( 32 1 amplifies a tiny error. In conse- 



quence, even when we have observed that pi is properly 
described by Eq. (23) for the case a = 0.10 [Fig [sl a)], 
the expression for the spike-count variance [Eq. ( |35[ )] is 
qualitatively good but approximate [note, however, that 
there is a significant reduction even for small values of a, 
and Eq. (35) has a moderate relative error for a = 0.1]. 
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The linear spike-count variance reduction expressed by 



Eq. ( 35 1 is unbounded, which is obviously unreasonable 



(it enables negative values for the spike-count variance). 
Higher-order terms should oppose this linear reduction, 
as can be observed in Fig. [6] for the case a — 0.10. 

Since the spike-count variance reduction for a given 
firing frequency arises from the interplay between the 
FPT statistics and the presence of correlations, it would 
be interesting to disentangle to what extent each effect 
contributes to the observed reduction. To analyze this 
question wc shuffled the spike train, which maintains the 
first-order statistics (FPT statistics) destroying ISI corre- 
lations. The spike-count variance reduction for the shuf- 
fled spike train is shown in Fig. [6] as semifilled symbols. 
The theoretical analog corresponds to (An|,) given ex- 



clusively by Eq. ( 33 1 , since the sum of the SCC values is 
[1 + 2^1^-^pk = 1 instead of Eq. ( 32 )] , and it is shown 



as dashed lines in Fig. [6] As expected, theoretical and 
numerical results are in good accordance. The behav- 
ior observed for each contribution is reasonable. At low 
firing frequencies, within each ISI the exponential evo- 
lution of the adaptation process has decayed, meaning 
that each — a and correlations disappear. Corre- 
spondingly, the spike-count variance reduction is given 
exclusively by the FPT statistics (this case corresponds 
to Td — in US]), and we denote this regime as an 
intrinsic reduction in Fig. |6] At large firing frequencies, 
the adaptation process within a single ISI is essentially 
constant (i.e., oo). In this case, the collection of 

ISIs satisfies a quasistatic approximation |47| . and the 
spike-count variance of the shuffled spike train is indis- 
tinguishable from the case of no adaptation. However, in 
this hmit, correlations decay very slowly [(^q (1/ra) — > 1; 
see Fig. [Sjjc)], and even when each SCC is small because 
p\ is, the sum of correlations accumulates over many lags, 
building up a finite nonvanishing value. In Fig. [6j we de- 
note this range as a reduction due to correlations ( "corrs" 
in the figure). 



Equation ( 35 ) shows that the spike-count variance re- 



duction is independent of the firing frequency, for the 
perfect IF neuron model. As shown in the previous para- 
graph, in the limit of low as well as large firing frequen- 
cies, the mechanisms that account for the reduction are 
different. The fact that both effects infiuence the spike- 
count variance to the same extent, and furthermore, that 
these mechanisms exactly compensate each other at in- 
termediate firing frequencies is intriguing. Obviously, 
these results hold for this particular IF model. It would 
be interesting to analyze the contributions in other IF 
models; we expect analogous spike-count variance reduc- 
tions and the same limit behaviors, but not an indepen- 
dence on firing frequencies (in general, in other models 
the spike-count variance for a = depends on the firing 
frequency) . 



IV. DISCUSSION AND CONCLUDING 
REMARKS 



In this work we have analyzed the onset of correlations 
for a general neuron model, where an external input as 
well as an internal spike-based adaptation current drive 
the membrane potential. The external current is com- 
posed by a static input and fast fluctuations. For this 
system, the dependence of the adaptation current on the 
past history, through the initial states of the adaptation 
process, facilitates the development of correlations be- 
tween successive ISIs [SI [TO] . In the regime of 
slight adaptation, we have shown that correlations share 
a general structure across different models. By means 
of a hidden Markov model, we have explicitly derived 
the dependence of the SCC on different properties of 
the FPT statistics corresponding to the underlying time- 
inhomogeneous stochastic process, Eqs. (23) and (26 1. 



In this regime, for one-dimensional models such as IF 
neuron models, the necessary properties are given |431 - 
I46j . For any other (high-dimensional) model, whenever 
a slight exponential time-dependent current smoothly re- 
shapes the FPT statistics in comparison to the unper- 
turbed case (analogously to Fig. 1 in [13]), the expres- 
sions derived here apply. The geometric structure and 
exponential decay for the SCC is surprisingly simple and 
general for the scenario considered here. This kind of 
structure was first observed by Lindner and Schwalger 
for successive escapes of an overdamped Brownian parti- 
cle in a randomly modulated asymmetric double well, 
which can be modeled as transitions between discrete 
states [ini HI] . Posteriorly, these authors extended their 
results to a situation with multiple internal states |39j 
and, in particular, studied the case of negative corre- 
lations in neurons with inhibitory feedback (resembling 
current-based adaptation), with an appropriate scheme 
of transitions. Even when general and very promising, a 
quantitative evaluation of correlations in adapting neu- 
rons under this framework requires a procedure for the es- 
timation of transition rates, in a discrete version of adap- 
tation, obtained from dynamical models and/or experi- 
mental data. The results obtained here are in qualitative 
agreement to those obtained by Schwalger and Lindner 
[39], which implies that this estimation procedure could 
be an interesting topic to study. 

The development of negative correlations in successive 
ISIs in a spike train influences the encoding capabilities 
of a neural system [TH l2Tti24] . Here, we have analyzed 
the decrease in the asymptotic spike-count variance due 
to the intrinsic variability reduction of the FPT statis- 



tics and the presence of correlations [Eqs. (31|-(34)] in 
comparison to the case without adaptation. In particu- 
lar, for the perfect IF neuron model the decrease is ex- 
tremely simple and does not depend on the firing fre- 
quency, Eq. (35). This analysis is theoretically impor- 
tant, but it should be put in the proper context. The 
spike-count variance reduction given by Eq. pl[ ) is valid 
for the asymptotic limit, which can be unfeasible in real 
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neurons, especially at low firing rates. For systems oper- 
ating with finite temporal windows, the framework pre- 
sented here should be extended by using the complete for- 
malism derived by van Vreeswijk in '42], and finely used 
by Farkhooi et al. to analyze a population scheme |24j . 
In this case, the spike-count variance will be a function 
of the length of the temporal window used to compute 
the statistics, and obviously, the results presented here 
should agree in the limit of large windows. This behav- 
ior was outlined by Chacron et al. in a different version 
of adaptation (see results of the model without slow noise 
in Fig. 4 of [33]). That work also highlighted a possible 
explanation to the interesting finding made by Ratnam 
and Nelson [32] , where the spike-count variance exhibits a 
minimum as a function of the length of the counting win- 
dow (a possible behaviorally important phenomenon) . In 
their work, Chacron et al. demonstrated that a slow ex- 
ternal noise leaves relatively intact the short-range cor- 
relations, while destroying negative correlations at large 
lags, giving rise to small and slowly decaying positive 
correlations which dominate the asymptotic regime. On 



the other hand, across different neural systems it has 
been observed that the only significant SCC corresponds 
to the lag 1 [ini [23], which reinforces the idea that a 
slow external noise would be an important part of the 
incoming signal, in addition to fast fluctuations. Theo- 
retically, the analysis of the spike-count variance for IF 
neuron models driven by slow fluctuations were success- 
fully carried out via a quasistatic approximation ^34. ,36] . 
This suggests that the hidden Markov model used here 
to model negative correlations could be extended with 
a similar quasistatic approximation in order to include 
slow fluctuations in the analysis. 
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